swedes <- read.dta(file="swedes_wide.dta")

#POLYCORIC CORRELATIONS
library(polycor)

swedes2 <- swedes[swedes$attityd381_na==0,]
MZ <- rep(NA,1000)
DZ <- rep(NA,1000)
DZO <- rep(NA,1000)
DIFF <- rep(NA,1000)
for(i in 1:1000){
swedes3 <- swedes2[sample(dim(swedes2)[1],replace=T),]
MZ[i] <- polychor(swedes3[swedes3$bestzyg1==1,3],swedes3[swedes3$bestzyg2==1,21])
DZ[i] <- polychor(swedes3[swedes3$bestzyg1==2,3],swedes3[swedes3$bestzyg2==2,21])
DIFF[i] <- MZ[i] - DZ[i]}

c(mean(MZ),quantile(MZ,c(0.025,0.975)))
c(mean(DZ),quantile(DZ,c(0.025,0.975)))
c(mean(DZO),quantile(DZO,c(0.025,0.975)))
c(mean(DIFF),quantile(DIFF,c(0.025,0.975)))
CORR[i,14] <- sum(DIFF <= 0)/1000
table(swedes2$bestzyg1)[1]
table(swedes2$bestzyg1)[2]

